Functional Analysis
Setup
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(clusterProfiler)##
## clusterProfiler v4.10.1 For help: https://yulab-smu.top/biomedical-knowledge-mining-book/
##
## If you use clusterProfiler in published research, please cite:
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021, 2(3):100141
##
## Attaching package: 'clusterProfiler'
##
## The following object is masked from 'package:purrr':
##
## simplify
##
## The following object is masked from 'package:stats':
##
## filterlibrary(org.Hs.eg.db)## Loading required package: AnnotationDbi
## Loading required package: stats4
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
##
## The following objects are masked from 'package:lubridate':
##
## intersect, setdiff, union
##
## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union
##
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
##
## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which.max, which.min
##
## Loading required package: Biobase
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
##
## Loading required package: IRanges
## Loading required package: S4Vectors
##
## Attaching package: 'S4Vectors'
##
## The following object is masked from 'package:clusterProfiler':
##
## rename
##
## The following objects are masked from 'package:lubridate':
##
## second, second<-
##
## The following objects are masked from 'package:dplyr':
##
## first, rename
##
## The following object is masked from 'package:tidyr':
##
## expand
##
## The following object is masked from 'package:utils':
##
## findMatches
##
## The following objects are masked from 'package:base':
##
## expand.grid, I, unname
##
##
## Attaching package: 'IRanges'
##
## The following object is masked from 'package:clusterProfiler':
##
## slice
##
## The following object is masked from 'package:lubridate':
##
## %within%
##
## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice
##
## The following object is masked from 'package:purrr':
##
## reduce
##
## The following object is masked from 'package:grDevices':
##
## windows
##
##
## Attaching package: 'AnnotationDbi'
##
## The following object is masked from 'package:clusterProfiler':
##
## select
##
## The following object is masked from 'package:dplyr':
##
## selectlibrary(enrichplot)
library(pathview)## ##############################################################################
## Pathview is an open source software package distributed under GNU General
## Public License version 3 (GPLv3). Details of GPLv3 is available at
## http://www.gnu.org/licenses/gpl-3.0.html. Particullary, users are required to
## formally cite the original Pathview paper (not just mention it) in publications
## or products. For details, do citation("pathview") within R.
##
## The pathview downloads and uses KEGG data. Non-academic uses may require a KEGG
## license agreement (details at http://www.kegg.jp/kegg/legal.html).
## ##############################################################################library(DOSE)## DOSE v3.28.2 For help: https://yulab-smu.top/biomedical-knowledge-mining-book/
##
## If you use DOSE in published research, please cite:
## Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology Semantic and Enrichment analysis. Bioinformatics 2015, 31(4):608-609knitr::opts_chunk$set(fig.width = 10, dpi = 300)
dir.create(path = "./results/FA", showWarnings = FALSE)Results
TLR10 light vs. dark 3h
Gene ontology
out.dir <- "./results/FA/TLR10_light_3h_vs_TLR10_dark_3h/"
dir.create(out.dir)
res.TLR3h <- read.csv("results/DESeq/TLR10_light_3h_vs_TLR10_dark_3h_unfiltered.csv") %>%
as_tibble() %>%
arrange(desc(log2FoldChange))
gene_list.TLR3h <- res.TLR3h$log2FoldChange
names(gene_list.TLR3h) <- res.TLR3h$ENSG
set.seed(5892)
gseGO.TLR3h <- gseGO(gene_list.TLR3h,
OrgDb = org.Hs.eg.db,
keyType = "ENSEMBL",
maxGSSize = 500,
minGSSize = 15,
seed = TRUE)## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...gseGO.TLR3h <- setReadable(gseGO.TLR3h, "org.Hs.eg.db", "ENSEMBL")
gseGO.TLR3h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism Homo sapiens
## #...@setType BP
## #...@keytype ENSEMBL
## #...@geneList Named num [1:32287] 1.61 1.54 1.4 1.37 1.31 ...
## - attr(*, "names")= chr [1:32287] "ENSG00000204388" "ENSG00000204389" "ENSG00000051108" "ENSG00000175197" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...49 enriched terms found
## 'data.frame': 49 obs. of 11 variables:
## $ ID : chr "GO:0006457" "GO:0051604" "GO:0002181" "GO:0061077" ...
## $ Description : chr "protein folding" "protein maturation" "cytoplasmic translation" "chaperone-mediated protein folding" ...
## $ setSize : int 200 495 154 71 134 167 156 41 121 98 ...
## $ enrichmentScore: num 0.676 0.525 0.657 0.769 0.671 ...
## $ NES : num 2.01 1.66 1.9 2.02 1.92 ...
## $ pvalue : num 2.00e-09 8.10e-08 5.01e-07 1.27e-06 2.65e-06 ...
## $ p.adjust : num 9.65e-06 1.95e-04 8.05e-04 1.52e-03 2.55e-03 ...
## $ qvalue : num 8.33e-06 1.68e-04 6.95e-04 1.32e-03 2.20e-03 ...
## $ rank : num 3165 4282 3050 2216 3646 ...
## $ leading_edge : chr "tags=44%, list=10%, signal=40%" "tags=35%, list=13%, signal=31%" "tags=60%, list=9%, signal=55%" "tags=44%, list=7%, signal=41%" ...
## $ core_enrichment: chr "HSPA1B/HSPA1A/HSPA1L/DNAJB1/HSPB1/HSPH1/PPID/DNAJA1/GRPEL2/CHORDC1/ENGASE/DNAJB2/HSPA5/ARL2/AIP/LMAN2L/ST13/CDC"| __truncated__ "HSPA1B/HSPA1A/HSPA1L/DNAJB1/CHAC1/HSPB1/HSPH1/SIRT4/PPID/DNAJA1/GRPEL2/CHORDC1/ENGASE/DNAJB2/HSPA5/ARL2/AIP/NOL"| __truncated__ "EIF4A2/RPL28/RPS19/AARS1/UBA52/EIF3H/RPL30/RPL13A/RPS12/RPS4X/RPL31/RPLP1/RPS28/RPS15A/RPL37A/RPS9/RPL14/RPL13/"| __truncated__ "HSPA1B/HSPA1A/HSPA1L/DNAJB1/HSPB1/HSPH1/PPID/CHORDC1/DNAJB2/HSPA5/ST13/SDF2L1/SDF2/FKBP2/PFDN5/PPIB/FKBP4/PFDN4"| __truncated__ ...
## #...Citation
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
## clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
## The Innovation. 2021, 2(3):100141write.csv(data.frame(gseGO.TLR3h), file = paste0(out.dir, "gsea_go.csv"), row.names = F)
TLR3h.go.plot <- nrow(gseGO.TLR3h) > 0dotplot(gseGO.TLR3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
cnetplot(gseGO.TLR3h, showCategory = 5, color.params = list(foldChange = gene_list.TLR3h),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.## Warning: ggrepel: 5 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
gseGO.TLR3h.pwts <- pairwise_termsim(gseGO.TLR3h)
emapplot(gseGO.TLR3h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
KEGG (Kyoto Encyclopedia of Genes and Genomes)
ids <- bitr(names(gene_list.TLR3h), fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)## 'select()' returned 1:many mapping between keys and columnsids <- ids[!duplicated(ids$ENSEMBL),]
gene_list.TLR3h.entrez <- gene_list.TLR3h[names(gene_list.TLR3h) %in% ids$ENSEMBL]
names(gene_list.TLR3h.entrez) <- ids$ENTREZID
gene_list.TLR3h.entrez <- gene_list.TLR3h.entrez[!duplicated(names(gene_list.TLR3h.entrez))]set.seed(5923)
gseKEGG.TLR3h <- gseKEGG(gene_list.TLR3h.entrez,
organism = "hsa",
keyType = "ncbi-geneid",
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = T)## Reading KEGG annotation online: "https://rest.kegg.jp/link/hsa/pathway"...## Reading KEGG annotation online: "https://rest.kegg.jp/list/pathway/hsa"...## Reading KEGG annotation online: "https://rest.kegg.jp/conv/ncbi-geneid/hsa"...## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...gseKEGG.TLR3h <- setReadable(gseKEGG.TLR3h, "org.Hs.eg.db", "ENTREZID")
gseKEGG.TLR3h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism hsa
## #...@setType KEGG
## #...@keytype ENTREZID
## #...@geneList Named num [1:22414] 1.61 1.54 1.4 1.37 1.31 ...
## - attr(*, "names")= chr [1:22414] "3304" "3303" "9709" "1649" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...73 enriched terms found
## 'data.frame': 73 obs. of 11 variables:
## $ ID : chr "hsa03010" "hsa05171" "hsa04141" "hsa05020" ...
## $ Description : chr "Ribosome" "Coronavirus disease - COVID-19" "Protein processing in endoplasmic reticulum" "Prion disease" ...
## $ setSize : int 129 203 164 255 57 131 255 149 54 292 ...
## $ enrichmentScore: num 0.786 0.597 0.63 0.567 0.791 ...
## $ NES : num 2.49 1.98 2.04 1.93 2.22 ...
## $ pvalue : num 8.61e-17 3.21e-09 6.93e-09 6.23e-09 5.65e-08 ...
## $ p.adjust : num 2.87e-14 5.37e-07 5.79e-07 5.79e-07 3.77e-06 ...
## $ qvalue : num 1.87e-14 3.48e-07 3.76e-07 3.76e-07 2.45e-06 ...
## $ rank : num 1880 2592 2319 2411 2044 ...
## $ leading_edge : chr "tags=62%, list=8%, signal=57%" "tags=41%, list=12%, signal=37%" "tags=36%, list=10%, signal=32%" "tags=34%, list=11%, signal=30%" ...
## $ core_enrichment: chr "RPL28/RPS19/UBA52/RPL30/MRPL34/RPL13A/RPS12/RPS4X/RPL31/RPLP1/RPS28/RPS15A/RPL37A/RPS9/RPL14/RPL13/RPL37/RPS11/"| __truncated__ "RPL28/RPS19/UBA52/RPL30/RPL13A/RPS12/RPS4X/RPL31/RPLP1/RPS28/RPS15A/RPL37A/RPS9/RPL14/RPL13/RPL37/RPS11/RSL24D1"| __truncated__ "HSPA1B/HSPA1A/HERPUD1/DDIT3/HSPA1L/DNAJB1/ATF4/HSPH1/XBP1/NFE2L2/DNAJA1/DNAJB2/HSPA5/TMEM258/SSR4/UBE2D1/LMAN2/"| __truncated__ "HSPA1B/HSPA1A/DDIT3/HSPA1L/ATF4/HSPA5/NDUFA3/COX5B/COX4I1/COX7C/UQCR11/UQCRB/NDUFS8/ATF2/ATP5F1E/NDUFA2/UQCRQ/C"| __truncated__ ...
## #...Citation
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
## clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
## The Innovation. 2021, 2(3):100141write.csv(data.frame(gseKEGG.TLR3h), file = paste0(out.dir, "gsea_kegg.csv"), row.names = F)
TLR3h.kegg.plot <- nrow(gseKEGG.TLR3h) > 0dotplot(gseKEGG.TLR3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
cnetplot(gseKEGG.TLR3h, showCategory = 5, color.params = list(foldChange = gene_list.TLR3h.entrez),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
gseKEGG.TLR3h.pwts <- pairwise_termsim(gseKEGG.TLR3h)
emapplot(gseKEGG.TLR3h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
pathview(gene.data = gene_list.TLR3h.entrez, pathway.id = gseKEGG.TLR3h@result$ID[1], species = "hsa")## 'select()' returned 1:1 mapping between keys and columns## Info: Working in directory C:/Users/Daniel Zimmermann/Desktop/AnS_DESeq## Info: Writing image file hsa03010.pathview.pngfilename <- paste0(gseKEGG.TLR3h@result$ID[1], ".pathview.png")
file.rename(filename, paste0(out.dir, filename))## [1] TRUEknitr::include_graphics(paste0(out.dir, filename))
Disease ontology
set.seed(5534)
gseDO.TLR3h <- gseDO(gene_list.TLR3h.entrez,
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = T)## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...gseDO.TLR3h <- setReadable(gseDO.TLR3h, "org.Hs.eg.db", "ENTREZID")
gseDO.TLR3h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism Homo sapiens
## #...@setType DO
## #...@keytype ENTREZID
## #...@geneList Named num [1:22414] 1.61 1.54 1.4 1.37 1.31 ...
## - attr(*, "names")= chr [1:22414] "3304" "3303" "9709" "1649" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...21 enriched terms found
## 'data.frame': 21 obs. of 11 variables:
## $ ID : chr "DOID:1470" "DOID:2841" "DOID:12306" "DOID:1176" ...
## $ Description : chr "major depressive disorder" "asthma" "vitiligo" "bronchial disease" ...
## $ setSize : int 37 255 70 264 19 17 111 24 158 300 ...
## $ enrichmentScore: num 0.813 0.512 0.693 0.499 0.902 ...
## $ NES : num 2.11 1.73 2.01 1.7 2.08 ...
## $ pvalue : num 3.82e-06 4.72e-06 1.01e-05 2.23e-05 4.03e-05 ...
## $ p.adjust : num 0.00226 0.00226 0.00323 0.00534 0.00594 ...
## $ qvalue : num 0.00205 0.00205 0.00293 0.00485 0.00538 ...
## $ rank : num 1299 1977 1359 1977 583 ...
## $ leading_edge : chr "tags=14%, list=6%, signal=13%" "tags=13%, list=9%, signal=12%" "tags=23%, list=6%, signal=22%" "tags=12%, list=9%, signal=12%" ...
## $ core_enrichment: chr "HSPA1A/HSPA1L/COMT/PTEN/EIF4B" "HSPA1B/HSPA1A/HSPB1/NFE2L2/LTA4H/UCN/PRMT2/CAT/GSTP1/SOD1/GSTM3/SELPLG/IFNGR1/MIF/ORMDL1/HLA-C/ARG2/NQO1/ANG/LG"| __truncated__ "HSPA1A/XBP1/NFE2L2/CAT/COMT/GSTP1/TSLP/HLA-A/PSMB8/GPX1/HLA-C/PSMB9/MSRB2/HLA-B/FAS/PTGS2" "HSPA1B/HSPA1A/HSPB1/NFE2L2/LTA4H/UCN/PRMT2/CAT/GSTP1/SOD1/GSTM3/SELPLG/IFNGR1/MIF/ORMDL1/HLA-C/ARG2/NQO1/ANG/LG"| __truncated__ ...
## #...Citation
## Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an
## R/Bioconductor package for Disease Ontology Semantic and Enrichment
## analysis. Bioinformatics 2015, 31(4):608-609write.csv(data.frame(gseDO.TLR3h), file = paste0(out.dir, "gsea_do.csv"), row.names = F)
TLR3h.DO.plot <- nrow(gseDO.TLR3h) > 0dotplot(gseDO.TLR3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
cnetplot(gseDO.TLR3h, showCategory = 5, layout = "dh", color.params = list(foldChange = gene_list.TLR3h.entrez),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
gseDO.TLR3h.pwts <- pairwise_termsim(gseDO.TLR3h)
emapplot(gseDO.TLR3h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
TLR10 light vs. dark 16h
Gene ontology
out.dir <- "./results/FA/TLR10_light_16h_vs_TLR10_dark_16h/"
dir.create(out.dir)
res.TLR16h <- read.csv("results/DESeq/TLR10_light_16h_vs_TLR10_dark_16h_unfiltered.csv") %>%
as_tibble() %>%
arrange(desc(log2FoldChange))
gene_list.TLR16h <- res.TLR16h$log2FoldChange
names(gene_list.TLR16h) <- res.TLR16h$ENSG
set.seed(4735)
gseGO.TLR16h <- gseGO(gene_list.TLR16h,
OrgDb = org.Hs.eg.db,
keyType = "ENSEMBL",
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = TRUE)## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...gseGO.TLR16h <- setReadable(gseGO.TLR16h, "org.Hs.eg.db", "ENSEMBL")
gseGO.TLR16h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism Homo sapiens
## #...@setType BP
## #...@keytype ENSEMBL
## #...@geneList Named num [1:32287] 1.7 1.63 1.55 1.54 1.46 ...
## - attr(*, "names")= chr [1:32287] "ENSG00000140297" "ENSG00000108448" "ENSG00000234719" "ENSG00000087842" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...401 enriched terms found
## 'data.frame': 401 obs. of 11 variables:
## $ ID : chr "GO:0007059" "GO:0098813" "GO:0000280" "GO:0048285" ...
## $ Description : chr "chromosome segregation" "nuclear chromosome segregation" "nuclear division" "organelle fission" ...
## $ setSize : int 404 297 413 459 417 220 435 477 181 322 ...
## $ enrichmentScore: num -0.761 -0.78 -0.733 -0.712 -0.714 ...
## $ NES : num -2.32 -2.33 -2.25 -2.19 -2.19 ...
## $ pvalue : num 7.68e-35 2.55e-29 7.75e-29 1.62e-28 2.86e-26 ...
## $ p.adjust : num 3.70e-31 6.13e-26 1.24e-25 1.95e-25 2.76e-23 ...
## $ qvalue : num 3.15e-31 5.22e-26 1.06e-25 1.66e-25 2.35e-23 ...
## $ rank : num 2472 2572 2779 3019 3039 ...
## $ leading_edge : chr "tags=47%, list=8%, signal=44%" "tags=47%, list=8%, signal=44%" "tags=41%, list=9%, signal=38%" "tags=40%, list=9%, signal=37%" ...
## $ core_enrichment: chr "KIF2A/KMT5A/GOLGA8A/MND1/CEP63/ACTL6A/ARID1A/SMARCC2/MSH5/CHMP2B/DDB1/SLC25A5/HAUS2/DIS3L2/PUM2/NUP62/NR3C1/RB1"| __truncated__ "ACTR3/BAG6/KIF2A/KMT5A/MND1/ACTL6A/ARID1A/SMARCC2/MSH5/CHMP2B/DDB1/DIS3L2/NUP62/RB1/SMARCE1/CDK5RAP2/SEH1L/HNRN"| __truncated__ "MLH1/MAPRE1/RAD1/CDC42/EPS8/CHEK1/ACTR3/BAG6/KIF2A/KMT5A/MND1/MSH5/CHMP2B/DDB1/RAD51D/DIS3L2/NUP62/RB1/CKS2/REE"| __truncated__ "TESMIN/RAB11A/REEP3/TENT4A/RALBP1/GDAP1/ARHGEF10/MLH1/MAPRE1/RAD1/CDC42/EPS8/CORO1C/CHEK1/ACTR3/BAG6/MYO19/KIF2"| __truncated__ ...
## #...Citation
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
## clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
## The Innovation. 2021, 2(3):100141write.csv(data.frame(gseGO.TLR16h), file = paste0(out.dir, "gsea_go.csv"), row.names = F)
TLR16h.go.plot <- nrow(gseGO.TLR16h) > 0dotplot(gseGO.TLR16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
cnetplot(gseGO.TLR16h, showCategory = 5, color.params = list(foldChange = gene_list.TLR16h),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.## Warning: ggrepel: 3 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
gseGO.TLR16h.pwts <- pairwise_termsim(gseGO.TLR16h)
emapplot(gseGO.TLR16h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
KEGG
ids <- bitr(names(gene_list.TLR16h), fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)## 'select()' returned 1:many mapping between keys and columnsids <- ids[!duplicated(ids$ENSEMBL),]
gene_list.TLR16h.entrez <- gene_list.TLR16h[names(gene_list.TLR16h) %in% ids$ENSEMBL]
names(gene_list.TLR16h.entrez) <- ids$ENTREZID
gene_list.TLR16h.entrez <- gene_list.TLR16h.entrez[!duplicated(names(gene_list.TLR16h.entrez))]set.seed(5923)
gseKEGG.TLR16h <- gseKEGG(gene_list.TLR16h.entrez,
organism = "hsa",
keyType = "ncbi-geneid",
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = T)## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...gseKEGG.TLR16h <- setReadable(gseKEGG.TLR16h, "org.Hs.eg.db", "ENTREZID")
gseKEGG.TLR16h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism hsa
## #...@setType KEGG
## #...@keytype ENTREZID
## #...@geneList Named num [1:22414] 1.7 1.55 1.54 1.46 1.42 ...
## - attr(*, "names")= chr [1:22414] "9245" "729978" "8544" "2069" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...42 enriched terms found
## 'data.frame': 42 obs. of 11 variables:
## $ ID : chr "hsa04110" "hsa03010" "hsa04510" "hsa03030" ...
## $ Description : chr "Cell cycle" "Ribosome" "Focal adhesion" "DNA replication" ...
## $ setSize : int 155 129 195 36 41 54 173 212 106 38 ...
## $ enrichmentScore: num -0.753 0.664 -0.621 -0.857 -0.843 ...
## $ NES : num -2.27 2.12 -1.92 -2.1 -2.1 ...
## $ pvalue : num 4.04e-16 1.18e-09 6.94e-09 6.23e-08 6.54e-08 ...
## $ p.adjust : num 1.35e-13 1.98e-07 7.72e-07 3.12e-06 3.12e-06 ...
## $ qvalue : num 1.06e-13 1.56e-07 6.08e-07 2.46e-06 2.46e-06 ...
## $ rank : num 1963 6043 3758 2048 2101 ...
## $ leading_edge : chr "tags=47%, list=9%, signal=43%" "tags=71%, list=27%, signal=52%" "tags=42%, list=17%, signal=35%" "tags=78%, list=9%, signal=71%" ...
## $ core_enrichment: chr "RB1/CCND1/TGFB2/CDC25B/EP300/ANAPC15/PDS5B/STAG1/ANAPC5/CCNH/PTTG1/NIPBL/ATM/CDT1/CDKN2D/BUB3/CREBBP/WEE1/ATR/C"| __truncated__ "MRPL23/RPS27/MRPL34/MRPL2/MRPL10/RPS15A/RPS19/RSL24D1/MRPS17/RPL28/RPL39/MRPL9/RPS28/MRPL16/RPL13A/MRPL36/MRPL1"| __truncated__ "IGF1R/PDGFB/COL4A1/ITGB5/PAK2/ITGB3/LAMA5/COL6A6/EMP2/LAMA4/LAMB1/ITGA5/CAV2/GSK3B/PIP5K1A/ITGA10/MYL12A/TNXB/M"| __truncated__ "RPA3/RNASEH1/POLE3/RFC4/RNASEH2A/POLD3/POLD1/RFC1/POLA2/RFC5/PCNA/PRIM1/RPA1/PRIM2/POLE2/MCM6/RFC2/FEN1/MCM4/MC"| __truncated__ ...
## #...Citation
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
## clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
## The Innovation. 2021, 2(3):100141write.csv(data.frame(gseKEGG.TLR16h), file = paste0(out.dir, "gsea_kegg.csv"), row.names = F)
TLR16h.kegg.plot <- nrow(gseKEGG.TLR16h) > 0dotplot(gseKEGG.TLR16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
cnetplot(gseKEGG.TLR16h, showCategory = 5, color.params = list(foldChange = gene_list.TLR16h.entrez),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
gseKEGG.TLR16h.pwts <- pairwise_termsim(gseKEGG.TLR16h)
emapplot(gseKEGG.TLR16h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
pathview(gene.data = gene_list.TLR16h.entrez, pathway.id = gseKEGG.TLR16h@result$ID[1], species = "hsa")## 'select()' returned 1:1 mapping between keys and columns## Info: Working in directory C:/Users/Daniel Zimmermann/Desktop/AnS_DESeq## Info: Writing image file hsa04110.pathview.pngfilename <- paste0(gseKEGG.TLR16h@result$ID[1], ".pathview.png")
file.rename(filename, paste0(out.dir, filename))## [1] TRUEknitr::include_graphics(paste0(out.dir, filename))
Disease ontology
set.seed(6246)
gseDO.TLR16h <- gseDO(gene_list.TLR16h.entrez,
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = T)## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...gseDO.TLR16h <- setReadable(gseDO.TLR16h, "org.Hs.eg.db", "ENTREZID")
gseDO.TLR16h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism Homo sapiens
## #...@setType DO
## #...@keytype ENTREZID
## #...@geneList Named num [1:22414] 1.7 1.55 1.54 1.46 1.42 ...
## - attr(*, "names")= chr [1:22414] "9245" "729978" "8544" "2069" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...96 enriched terms found
## 'data.frame': 96 obs. of 11 variables:
## $ ID : chr "DOID:0080015" "DOID:5683" "DOID:2490" "DOID:1793" ...
## $ Description : chr "physical disorder" "hereditary breast ovarian cancer syndrome" "congenital nervous system abnormality" "pancreatic cancer" ...
## $ setSize : int 302 213 135 453 460 264 255 144 144 29 ...
## $ enrichmentScore: num -0.581 -0.617 -0.665 -0.493 -0.489 ...
## $ NES : num -1.87 -1.91 -1.97 -1.63 -1.62 ...
## $ pvalue : num 7.25e-10 5.94e-09 1.94e-08 3.13e-07 2.54e-07 ...
## $ p.adjust : num 6.96e-07 2.85e-06 6.19e-06 6.00e-05 6.00e-05 ...
## $ qvalue : num 5.29e-07 2.17e-06 4.71e-06 4.57e-05 4.57e-05 ...
## $ rank : num 3211 2757 3211 2639 2222 ...
## $ leading_edge : chr "tags=33%, list=14%, signal=29%" "tags=40%, list=12%, signal=36%" "tags=40%, list=14%, signal=34%" "tags=26%, list=12%, signal=24%" ...
## $ core_enrichment: chr "NEDD4L/NOG/LAMB1/HMGA2/ARFGEF2/GOLGB1/STAG2/FGFR1/VANGL1/NUAK2/MAD1L1/MEST/CEP290/NSD1/DAW1/MECP2/TWSG1/RPGRIP1"| __truncated__ "CCNT1/NOLC1/EMSY/NPM1/MLH1/CHD8/CHEK1/GTF2I/DHX9/SP1/FASN/PALB2/NCOA2/RAD51D/HNRNPA2B1/COL1A1/RB1/CCND1/ERCC5/R"| __truncated__ "NEDD4L/LAMB1/ARFGEF2/STAG2/NUAK2/MAD1L1/NSD1/MECP2/TWSG1/RPGRIP1L/PNKP/CHEK1/DAG1/MTRR/MAP1B/CDK5RAP2/AUTS2/REL"| __truncated__ "RALBP1/PRKD1/NRP2/FANCC/MLH1/IKBKB/TSC2/LIG4/CDC42/BRAF/CHEK1/MTRR/APAF1/ILK/MAPK9/DNMT3B/ANXA2/CYCS/CTNNA1/GDI"| __truncated__ ...
## #...Citation
## Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an
## R/Bioconductor package for Disease Ontology Semantic and Enrichment
## analysis. Bioinformatics 2015, 31(4):608-609write.csv(data.frame(gseDO.TLR16h), file = paste0(out.dir, "gsea_do.csv"), row.names = F)
TLR16h.DO.plot <- nrow(gseDO.TLR16h) > 0dotplot(gseDO.TLR16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
cnetplot(gseDO.TLR16h, showCategory = 5, color.params = list(foldChange = gene_list.TLR16h.entrez),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.## Warning: ggrepel: 13 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
gseDO.TLR16h.pwts <- pairwise_termsim(gseDO.TLR16h)
emapplot(gseDO.TLR16h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
TLR10 light vs. wildtype light 3h
Gene ontology
out.dir <- "./results/FA/TLR10_light_3h_vs_wildtype_light_3h/"
dir.create(out.dir)
res.TLRlightWT3h <- read.csv("results/DESeq/TLR10_light_3h_vs_wildtype_light_3h_unfiltered.csv") %>%
as_tibble() %>%
arrange(desc(log2FoldChange))
gene_list.TLRlightWT3h <- res.TLRlightWT3h$log2FoldChange
names(gene_list.TLRlightWT3h) <- res.TLRlightWT3h$ENSG
set.seed(4735)
gseGO.TLRlightWT3h <- gseGO(gene_list.TLRlightWT3h,
OrgDb = org.Hs.eg.db,
keyType = "ENSEMBL",
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = TRUE)## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...gseGO.TLRlightWT3h <- setReadable(gseGO.TLRlightWT3h, "org.Hs.eg.db", "ENSEMBL")
gseGO.TLRlightWT3h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism Homo sapiens
## #...@setType BP
## #...@keytype ENSEMBL
## #...@geneList Named num [1:32287] 9.67 3.47 2.11 1.99 1.96 ...
## - attr(*, "names")= chr [1:32287] "ENSG00000174123" "ENSG00000164185" "ENSG00000123453" "ENSG00000225756" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...12 enriched terms found
## 'data.frame': 12 obs. of 11 variables:
## $ ID : chr "GO:0030865" "GO:0002224" "GO:0030866" "GO:0032528" ...
## $ Description : chr "cortical cytoskeleton organization" "toll-like receptor signaling pathway" "cortical actin cytoskeleton organization" "microvillus organization" ...
## $ setSize : int 50 63 42 24 58 213 83 65 56 68 ...
## $ enrichmentScore: num -0.997 0.998 -0.998 -0.999 -0.996 ...
## $ NES : num -2.27 2.16 -2.22 -2.08 -2.29 ...
## $ pvalue : num 1.15e-06 1.63e-06 4.89e-06 5.39e-06 1.88e-05 ...
## $ p.adjust : num 0.00393 0.00393 0.00649 0.00649 0.01634 ...
## $ qvalue : num 0.0037 0.0037 0.00612 0.00612 0.0154 ...
## $ rank : num 2 7 2 10 2 2 2 6 6 6 ...
## $ leading_edge : chr "tags=4%, list=0%, signal=4%" "tags=3%, list=0%, signal=3%" "tags=5%, list=0%, signal=5%" "tags=12%, list=0%, signal=13%" ...
## $ core_enrichment: chr "PLEK/LCP1" "TLR10/CD40" "PLEK/LCP1" "GLDN/TNIK/STK26" ...
## #...Citation
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
## clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
## The Innovation. 2021, 2(3):100141write.csv(data.frame(gseGO.TLRlightWT3h), file = paste0(out.dir, "gsea_go.csv"), row.names = F)
TLRlightWT3h.go.plot <- nrow(gseGO.TLRlightWT3h) > 0dotplot(gseGO.TLRlightWT3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
cnetplot(gseGO.TLRlightWT3h, showCategory = 5, color.params = list(foldChange = gene_list.TLRlightWT3h),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
gseGO.TLRlightWT3h.pwts <- pairwise_termsim(gseGO.TLRlightWT3h)
emapplot(gseGO.TLRlightWT3h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
KEGG
ids <- bitr(names(gene_list.TLRlightWT3h), fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)## 'select()' returned 1:many mapping between keys and columnsids <- ids[!duplicated(ids$ENSEMBL),]
gene_list.TLRlightWT3h.entrez <- gene_list.TLRlightWT3h[names(gene_list.TLRlightWT3h) %in% ids$ENSEMBL]
names(gene_list.TLRlightWT3h.entrez) <- ids$ENTREZID
gene_list.TLRlightWT3h.entrez <- gene_list.TLRlightWT3h.entrez[!duplicated(names(gene_list.TLRlightWT3h.entrez))]set.seed(5923)
gseKEGG.TLRlightWT3h <- gseKEGG(gene_list.TLRlightWT3h.entrez,
organism = "hsa",
keyType = "ncbi-geneid",
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = T)## preparing geneSet collections...## GSEA analysis...## no term enriched under specific pvalueCutoff...gseKEGG.TLRlightWT3h <- setReadable(gseKEGG.TLRlightWT3h, "org.Hs.eg.db", "ENTREZID")
gseKEGG.TLRlightWT3h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism hsa
## #...@setType KEGG
## #...@keytype ENTREZID
## #...@geneList Named num [1:22414] 9.67 3.47 2.11 1.99 1.96 ...
## - attr(*, "names")= chr [1:22414] "81793" "133923" "1757" "138948" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...0 enriched terms found
## 'data.frame': 0 obs. of 8 variables:
## $ ID : chr
## $ Description : chr
## $ setSize : int
## $ enrichmentScore: num
## $ NES : num
## $ pvalue : num
## $ p.adjust : num
## $ qvalue : num
## #...Citation
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
## clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
## The Innovation. 2021, 2(3):100141write.csv(data.frame(gseKEGG.TLRlightWT3h), file = paste0(out.dir, "gsea_kegg.csv"), row.names = F)
TLRlightWT3h.kegg.plot <- nrow(gseKEGG.TLRlightWT3h) > 0dotplot(gseKEGG.TLRlightWT3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")cnetplot(gseKEGG.TLRlightWT3h, showCategory = 5, color.params = list(foldChange = gene_list.TLRlightWT3h.entrez),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)gseKEGG.TLRlightWT3h.pwts <- pairwise_termsim(gseKEGG.TLRlightWT3h)
emapplot(gseKEGG.TLRlightWT3h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")pathview(gene.data = gene_list.TLRlightWT3h.entrez, pathway.id = gseKEGG.TLRlightWT3h@result$ID[1], species = "hsa")
filename <- paste0(gseKEGG.TLRlightWT3h@result$ID[1], ".pathview.png")
file.rename(filename, paste0(out.dir, filename))
knitr::include_graphics(paste0(out.dir, filename))Disease ontology
set.seed(6246)
gseDO.TLRlightWT3h <- gseDO(gene_list.TLRlightWT3h.entrez,
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = T)## preparing geneSet collections...## GSEA analysis...## no term enriched under specific pvalueCutoff...gseDO.TLRlightWT3h <- setReadable(gseDO.TLRlightWT3h, "org.Hs.eg.db", "ENTREZID")
gseDO.TLRlightWT3h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism Homo sapiens
## #...@setType DO
## #...@keytype ENTREZID
## #...@geneList Named num [1:22414] 9.67 3.47 2.11 1.99 1.96 ...
## - attr(*, "names")= chr [1:22414] "81793" "133923" "1757" "138948" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...0 enriched terms found
## 'data.frame': 0 obs. of 8 variables:
## $ ID : chr
## $ Description : chr
## $ setSize : int
## $ enrichmentScore: num
## $ NES : num
## $ pvalue : num
## $ p.adjust : num
## $ qvalue : num
## #...Citation
## Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an
## R/Bioconductor package for Disease Ontology Semantic and Enrichment
## analysis. Bioinformatics 2015, 31(4):608-609write.csv(data.frame(gseDO.TLRlightWT3h), file = paste0(out.dir, "gsea_do.csv"), row.names = F)
TLRlightWT3h.DO.plot <- nrow(gseDO.TLRlightWT3h) > 0dotplot(gseDO.TLRlightWT3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")cnetplot(gseDO.TLRlightWT3h, showCategory = 5, color.params = list(foldChange = gene_list.TLRlightWT3h.entrez),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)gseDO.TLRlightWT3h.pwts <- pairwise_termsim(gseDO.TLRlightWT3h)
emapplot(gseDO.TLRlightWT3h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")TLR10 light vs. wildtype light 16h
Gene ontology
out.dir <- "./results/FA/TLR10_light_16h_vs_wildtype_light_16h/"
dir.create(out.dir)
res.TLRlightWT16h <- read.csv("results/DESeq/TLR10_light_16h_vs_wildtype_light_16h_unfiltered.csv") %>%
as_tibble() %>%
arrange(desc(log2FoldChange))
gene_list.TLRlightWT16h <- res.TLRlightWT16h$log2FoldChange
names(gene_list.TLRlightWT16h) <- res.TLRlightWT16h$ENSG
set.seed(6548)
gseGO.TLRlightWT16h <- gseGO(gene_list.TLRlightWT16h,
OrgDb = org.Hs.eg.db,
keyType = "ENSEMBL",
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = TRUE)## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...gseGO.TLRlightWT16h <- setReadable(gseGO.TLRlightWT16h, "org.Hs.eg.db", "ENSEMBL")
gseGO.TLRlightWT16h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism Homo sapiens
## #...@setType BP
## #...@keytype ENSEMBL
## #...@geneList Named num [1:32287] 6.72 2.57 2.19 1.92 1.57 ...
## - attr(*, "names")= chr [1:32287] "ENSG00000174123" "ENSG00000225756" "ENSG00000139292" "ENSG00000123453" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...3 enriched terms found
## 'data.frame': 3 obs. of 11 variables:
## $ ID : chr "GO:2000052" "GO:2000095" "GO:0072111"
## $ Description : chr "positive regulation of non-canonical Wnt signaling pathway" "regulation of Wnt signaling pathway, planar cell polarity pathway" "cell proliferation involved in kidney development"
## $ setSize : int 15 16 21
## $ enrichmentScore: num -0.986 -0.986 -0.987
## $ NES : num -1.87 -1.91 -2.01
## $ pvalue : num 1.66e-06 4.25e-06 1.12e-05
## $ p.adjust : num 0.00798 0.01024 0.01792
## $ qvalue : num 0.00678 0.0087 0.01523
## $ rank : num 2 148 2
## $ leading_edge : chr "tags=13%, list=0%, signal=13%" "tags=12%, list=0%, signal=12%" "tags=10%, list=0%, signal=10%"
## $ core_enrichment: chr "WNT5A/GPC3" "SFRP2/GPC3" "ITGB3/GPC3"
## #...Citation
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
## clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
## The Innovation. 2021, 2(3):100141write.csv(data.frame(gseGO.TLRlightWT16h), file = paste0(out.dir, "gsea_go.csv"), row.names = F)
TLRlightWT16h.go.plot <- nrow(gseGO.TLRlightWT16h) > 0dotplot(gseGO.TLRlightWT16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
cnetplot(gseGO.TLRlightWT16h, showCategory = 5, color.params = list(foldChange = gene_list.TLRlightWT16h),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
gseGO.TLRlightWT16h.pwts <- pairwise_termsim(gseGO.TLRlightWT16h)
emapplot(gseGO.TLRlightWT16h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
KEGG
ids <- bitr(names(gene_list.TLRlightWT16h), fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)## 'select()' returned 1:many mapping between keys and columnsids <- ids[!duplicated(ids$ENSEMBL),]
gene_list.TLRlightWT16h.entrez <- gene_list.TLRlightWT16h[names(gene_list.TLRlightWT16h) %in% ids$ENSEMBL]
names(gene_list.TLRlightWT16h.entrez) <- ids$ENTREZID
gene_list.TLRlightWT16h.entrez <- gene_list.TLRlightWT16h.entrez[!duplicated(names(gene_list.TLRlightWT16h.entrez))]set.seed(5923)
gseKEGG.TLRlightWT16h <- gseKEGG(gene_list.TLRlightWT16h.entrez,
organism = "hsa",
keyType = "ncbi-geneid",
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = T)## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...gseKEGG.TLRlightWT16h <- setReadable(gseKEGG.TLRlightWT16h, "org.Hs.eg.db", "ENTREZID")
gseKEGG.TLRlightWT16h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism hsa
## #...@setType KEGG
## #...@keytype ENTREZID
## #...@geneList Named num [1:22414] 6.72 2.57 2.19 1.92 1.57 ...
## - attr(*, "names")= chr [1:22414] "81793" "138948" "8549" "1757" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...1 enriched terms found
## 'data.frame': 1 obs. of 11 variables:
## $ ID : chr "hsa03010"
## $ Description : chr "Ribosome"
## $ setSize : int 129
## $ enrichmentScore: num 0.746
## $ NES : num 2.25
## $ pvalue : num 9.41e-06
## $ p.adjust : num 0.00314
## $ qvalue : num 0.00276
## $ rank : int 3268
## $ leading_edge : chr "tags=69%, list=15%, signal=59%"
## $ core_enrichment: chr "MRPS17/RPS19/MRPL34/RPS15A/RPS29/RPLP1/RPL28/UBA52/RPL13A/RPS27L/RPL13/MRPS21/RPL37A/MRPL36/RPL27A/RPLP2/MRPL4/"| __truncated__
## #...Citation
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
## clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
## The Innovation. 2021, 2(3):100141write.csv(data.frame(gseKEGG.TLRlightWT16h), file = paste0(out.dir, "gsea_kegg.csv"), row.names = F)
TLRlightWT16h.kegg.plot <- nrow(gseKEGG.TLRlightWT16h) > 0dotplot(gseKEGG.TLRlightWT16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
cnetplot(gseKEGG.TLRlightWT16h, showCategory = 5, color.params = list(foldChange = gene_list.TLRlightWT16h.entrez),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
gseKEGG.TLRlightWT16h.pwts <- pairwise_termsim(gseKEGG.TLRlightWT16h)
emapplot(gseKEGG.TLRlightWT16h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")
pathview(gene.data = gene_list.TLRlightWT16h.entrez, pathway.id = gseKEGG.TLRlightWT16h@result$ID[1], species = "hsa")## 'select()' returned 1:1 mapping between keys and columns## Info: Working in directory C:/Users/Daniel Zimmermann/Desktop/AnS_DESeq## Info: Writing image file hsa03010.pathview.pngfilename <- paste0(gseKEGG.TLRlightWT16h@result$ID[1], ".pathview.png")
file.rename(filename, paste0(out.dir, filename))## [1] TRUEknitr::include_graphics(paste0(out.dir, filename))
Disease ontology
set.seed(6246)
gseDO.TLRlightWT16h <- gseDO(gene_list.TLRlightWT16h.entrez,
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = T)## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...gseDO.TLRlightWT16h <- setReadable(gseDO.TLRlightWT16h, "org.Hs.eg.db", "ENTREZID")
gseDO.TLRlightWT16h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism Homo sapiens
## #...@setType DO
## #...@keytype ENTREZID
## #...@geneList Named num [1:22414] 6.72 2.57 2.19 1.92 1.57 ...
## - attr(*, "names")= chr [1:22414] "81793" "138948" "8549" "1757" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...33 enriched terms found
## 'data.frame': 33 obs. of 11 variables:
## $ ID : chr "DOID:3113" "DOID:4481" "DOID:2021" "DOID:3594" ...
## $ Description : chr "papillary carcinoma" "allergic rhinitis" "placenta cancer" "choriocarcinoma" ...
## $ setSize : int 44 87 33 33 154 118 51 63 33 128 ...
## $ enrichmentScore: num -0.962 0.813 -0.973 -0.973 0.706 ...
## $ NES : num -2.41 2.36 -2.32 -2.32 2.2 ...
## $ pvalue : num 1.38e-05 1.14e-05 7.06e-06 7.06e-06 1.07e-05 ...
## $ p.adjust : num 0.00264 0.00264 0.00264 0.00264 0.00264 ...
## $ qvalue : num 0.00237 0.00237 0.00237 0.00237 0.00237 ...
## $ rank : num 2 5 2 2 5 5 333 2 2 5 ...
## $ leading_edge : chr "tags=5%, list=0%, signal=5%" "tags=2%, list=0%, signal=2%" "tags=6%, list=0%, signal=6%" "tags=6%, list=0%, signal=6%" ...
## $ core_enrichment: chr "MET/GPC3" "TLR10/CD40" "MTOR/GPC3" "MTOR/GPC3" ...
## #...Citation
## Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an
## R/Bioconductor package for Disease Ontology Semantic and Enrichment
## analysis. Bioinformatics 2015, 31(4):608-609write.csv(data.frame(gseDO.TLRlightWT16h), file = paste0(out.dir, "gsea_do.csv"), row.names = F)
TLRlightWT16h.DO.plot <- nrow(gseDO.TLRlightWT16h) > 0dotplot(gseDO.TLRlightWT16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
cnetplot(gseDO.TLRlightWT16h, showCategory = 5, color.params = list(foldChange = gene_list.TLRlightWT16h.entrez),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
gseDO.TLRlightWT16h.pwts <- pairwise_termsim(gseDO.TLRlightWT16h)
emapplot(gseDO.TLRlightWT16h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
TLR10 dark vs. wildtype light 3h
Gene ontology
out.dir <- "./results/FA/TLR10_dark_3h_vs_wildtype_light_3h/"
dir.create(out.dir)
res.TLRdarkWT3h <- read.csv("results/DESeq/TLR10_dark_3h_vs_wildtype_light_3h_unfiltered.csv") %>%
as_tibble() %>%
arrange(desc(log2FoldChange))
gene_list.TLRdarkWT3h <- res.TLRdarkWT3h$log2FoldChange
names(gene_list.TLRdarkWT3h) <- res.TLRdarkWT3h$ENSG
set.seed(4735)
gseGO.TLRdarkWT3h <- gseGO(gene_list.TLRdarkWT3h,
OrgDb = org.Hs.eg.db,
keyType = "ENSEMBL",
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = TRUE)## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...gseGO.TLRdarkWT3h <- setReadable(gseGO.TLRdarkWT3h, "org.Hs.eg.db", "ENSEMBL")
gseGO.TLRdarkWT3h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism Homo sapiens
## #...@setType BP
## #...@keytype ENSEMBL
## #...@geneList Named num [1:32287] 28.9 27.5 19.4 11.2 10.1 ...
## - attr(*, "names")= chr [1:32287] "ENSG00000157851" "ENSG00000132182" "ENSG00000165588" "ENSG00000130294" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...17 enriched terms found
## 'data.frame': 17 obs. of 11 variables:
## $ ID : chr "GO:0006986" "GO:0031667" "GO:0009991" "GO:0070059" ...
## $ Description : chr "response to unfolded protein" "response to nutrient levels" "response to extracellular stimulus" "intrinsic apoptotic signaling pathway in response to endoplasmic reticulum stress" ...
## $ setSize : int 134 437 464 61 200 202 365 495 296 170 ...
## $ enrichmentScore: num -0.764 -0.604 -0.599 -0.852 -0.702 ...
## $ NES : num -2 -1.78 -1.77 -2.01 -1.92 ...
## $ pvalue : num 3.37e-06 1.25e-06 2.42e-06 1.18e-05 1.05e-05 ...
## $ p.adjust : num 0.00542 0.00542 0.00542 0.01139 0.01139 ...
## $ qvalue : num 0.00526 0.00526 0.00526 0.01106 0.01106 ...
## $ rank : num 1924 3563 3563 767 3196 ...
## $ leading_edge : chr "tags=24%, list=6%, signal=23%" "tags=26%, list=11%, signal=24%" "tags=26%, list=11%, signal=24%" "tags=25%, list=2%, signal=24%" ...
## $ core_enrichment: chr "RHBDD1/DDRGK1/CREB3L4/ERN1/CDK5RAP3/HSPA6/PARP8/HSPA8/BFAR/OS9/CTH/EDEM2/ATF3/DNAJB2/CREBRF/RACK1/MANF/SERPINH1"| __truncated__ "SORL1/SPX/COMT/PRKAB2/TSSK6/MAP1LC3A/GABARAPL1/RXRA/LRRK2/VPS41/ULK1/MLST8/FES/APOE/PARP2/PPARA/PLD1/TSPO/ACE/L"| __truncated__ "SORL1/SPX/COMT/PRKAB2/TSSK6/MAP1LC3A/GABARAPL1/RXRA/NOCT/LRRK2/VPS41/ULK1/MLST8/FES/APOE/PARP2/PPARA/PLD1/TSPO/"| __truncated__ "HYOU1/CASP4/SYVN1/BCL2/BBC3/ERP29/CEBPB/TRIB3/PMAIP1/XBP1/ATF4/HSPA1A/HERPUD1/CHAC1/DDIT3" ...
## #...Citation
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
## clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
## The Innovation. 2021, 2(3):100141write.csv(data.frame(gseGO.TLRdarkWT3h), file = paste0(out.dir, "gsea_go.csv"), row.names = F)
TLRdarkWT3h.go.plot <- nrow(gseGO.TLRdarkWT3h) > 0dotplot(gseGO.TLRdarkWT3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
cnetplot(gseGO.TLRdarkWT3h, showCategory = 5, color.params = list(foldChange = gene_list.TLRdarkWT3h),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
gseGO.TLRdarkWT3h.pwts <- pairwise_termsim(gseGO.TLRdarkWT3h)
emapplot(gseGO.TLRdarkWT3h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
KEGG
ids <- bitr(names(gene_list.TLRdarkWT3h), fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)## 'select()' returned 1:many mapping between keys and columnsids <- ids[!duplicated(ids$ENSEMBL),]
gene_list.TLRdarkWT3h.entrez <- gene_list.TLRdarkWT3h[names(gene_list.TLRdarkWT3h) %in% ids$ENSEMBL]
names(gene_list.TLRdarkWT3h.entrez) <- ids$ENTREZID
gene_list.TLRdarkWT3h.entrez <- gene_list.TLRdarkWT3h.entrez[!duplicated(names(gene_list.TLRdarkWT3h.entrez))]set.seed(5923)
gseKEGG.TLRdarkWT3h <- gseKEGG(gene_list.TLRdarkWT3h.entrez,
organism = "hsa",
keyType = "ncbi-geneid",
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = T)## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...gseKEGG.TLRdarkWT3h <- setReadable(gseKEGG.TLRdarkWT3h, "org.Hs.eg.db", "ENTREZID")
gseKEGG.TLRdarkWT3h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism hsa
## #...@setType KEGG
## #...@keytype ENTREZID
## #...@geneList Named num [1:22414] 28.9 27.5 19.4 11.2 10.1 ...
## - attr(*, "names")= chr [1:22414] "56896" "23225" "5015" "547" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...3 enriched terms found
## 'data.frame': 3 obs. of 11 variables:
## $ ID : chr "hsa03010" "hsa04141" "hsa04612"
## $ Description : chr "Ribosome" "Protein processing in endoplasmic reticulum" "Antigen processing and presentation"
## $ setSize : int 129 164 57
## $ enrichmentScore: num -0.735 -0.696 -0.812
## $ NES : num -2.05 -2.01 -2.02
## $ pvalue : num 6.10e-06 3.56e-06 1.33e-04
## $ p.adjust : num 0.00102 0.00102 0.01475
## $ qvalue : num 0.000979 0.000979 0.014182
## $ rank : num 2463 2245 1726
## $ leading_edge : chr "tags=59%, list=11%, signal=53%" "tags=29%, list=10%, signal=27%" "tags=32%, list=8%, signal=29%"
## $ core_enrichment: chr "RPS15A/RPL36AL/RPL41/MRPS15/RPL24/RPL35/RPL39/RPL23A/MRPL10/RPS5/RPL10/RPL6/RPLP0/RPS27A/RPL15/RPS7/RPL22/RPL18"| __truncated__ "SSR2/DNAJC3/DNAJC1/TRAF2/UBXN6/TUSC3/LMAN1/ERO1B/MAN1B1/TXNDC5/PDIA4/BAG2/SIL1/ERN1/HSPA6/PRKCSH/CRYAB/DNAJB11/"| __truncated__ "NFYA/HLA-E/HLA-F/HSPA6/HSPA8/HLA-B/NFYB/B2M/CTSL/HLA-C/PSME1/HLA-A/CREB1/TAP1/HSPA5/HSPA1L/HSPA1A/HSPA1B"
## #...Citation
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
## clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
## The Innovation. 2021, 2(3):100141write.csv(data.frame(gseKEGG.TLRdarkWT3h), file = paste0(out.dir, "gsea_kegg.csv"), row.names = F)
TLRdarkWT3h.kegg.plot <- nrow(gseKEGG.TLRdarkWT3h) > 0dotplot(gseKEGG.TLRdarkWT3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
cnetplot(gseKEGG.TLRdarkWT3h, showCategory = 5, color.params = list(foldChange = gene_list.TLRdarkWT3h.entrez),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
gseKEGG.TLRdarkWT3h.pwts <- pairwise_termsim(gseKEGG.TLRdarkWT3h)
emapplot(gseKEGG.TLRdarkWT3h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
pathview(gene.data = gene_list.TLRdarkWT3h.entrez, pathway.id = gseKEGG.TLRdarkWT3h@result$ID[1], species = "hsa")## 'select()' returned 1:1 mapping between keys and columns## Info: Working in directory C:/Users/Daniel Zimmermann/Desktop/AnS_DESeq## Info: Writing image file hsa03010.pathview.pngfilename <- paste0(gseKEGG.TLRdarkWT3h@result$ID[1], ".pathview.png")
file.rename(filename, paste0(out.dir, filename))## [1] TRUEknitr::include_graphics(paste0(out.dir, filename))
Disease ontology
set.seed(6246)
gseDO.TLRdarkWT3h <- gseDO(gene_list.TLRdarkWT3h.entrez,
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = T)## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...gseDO.TLRdarkWT3h <- setReadable(gseDO.TLRdarkWT3h, "org.Hs.eg.db", "ENTREZID")
gseDO.TLRdarkWT3h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism Homo sapiens
## #...@setType DO
## #...@keytype ENTREZID
## #...@geneList Named num [1:22414] 28.9 27.5 19.4 11.2 10.1 ...
## - attr(*, "names")= chr [1:22414] "56896" "23225" "5015" "547" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...8 enriched terms found
## 'data.frame': 8 obs. of 11 variables:
## $ ID : chr "DOID:11949" "DOID:12306" "DOID:3021" "DOID:10825" ...
## $ Description : chr "Creutzfeldt-Jakob disease" "vitiligo" "acute kidney failure" "essential hypertension" ...
## $ setSize : int 23 70 99 111 289 459 79 470
## $ enrichmentScore: num -0.915 -0.783 -0.725 -0.711 -0.554 ...
## $ NES : num -2.01 -2 -1.93 -1.92 -1.68 ...
## $ pvalue : num 1.50e-04 1.93e-04 8.59e-05 1.31e-04 2.01e-04 ...
## $ p.adjust : num 0.0321 0.0321 0.0321 0.0321 0.0321 ...
## $ qvalue : num 0.0285 0.0285 0.0285 0.0285 0.0285 ...
## $ rank : num 243 2235 1728 2636 1861 ...
## $ leading_edge : chr "tags=13%, list=1%, signal=13%" "tags=24%, list=10%, signal=22%" "tags=22%, list=8%, signal=21%" "tags=23%, list=12%, signal=20%" ...
## $ core_enrichment: chr "ABCB1/SNCA/PTGS2" "HGF/IL6/FAS/TXNDC5/GSTP1/PSMB8/HLA-B/PSMB9/HLA-C/CAT/HLA-A/TAP1/TSLP/NFE2L2/XBP1/PTGS2/HSPA1A" "PPARG/GFER/LGALS3/FIS1/HMOX1/HSPA8/B2M/HYOU1/CLU/CASP1/SLC22A1/PTAFR/TNFRSF1B/CDKN1B/NQO1/BCL2/IL6R/NFE2L2/ARG2"| __truncated__ "FMO3/ADA/MMP2/GSTM3/ERAP1/TGFB3/HGF/IL6/ATP5PF/CST3/PPARG/HMOX1/MTHFR/KYNU/CTH/PSMB9/GDF15/ALDH2/CAT/CLU/RGS5/M"| __truncated__ ...
## #...Citation
## Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an
## R/Bioconductor package for Disease Ontology Semantic and Enrichment
## analysis. Bioinformatics 2015, 31(4):608-609write.csv(data.frame(gseDO.TLRdarkWT3h), file = paste0(out.dir, "gsea_do.csv"), row.names = F)
TLRdarkWT3h.DO.plot <- nrow(gseDO.TLRdarkWT3h) > 0dotplot(gseDO.TLRdarkWT3h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
cnetplot(gseDO.TLRdarkWT3h, showCategory = 5, color.params = list(foldChange = gene_list.TLRdarkWT3h.entrez),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
gseDO.TLRdarkWT3h.pwts <- pairwise_termsim(gseDO.TLRdarkWT3h)
emapplot(gseDO.TLRdarkWT3h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
TLR10 dark vs. wildtype light 16h
Gene ontology
out.dir <- "./results/FA/TLR10_dark_16h_vs_wildtype_light_16h/"
dir.create(out.dir)
res.TLRdarkWT16h <- read.csv("results/DESeq/TLR10_dark_16h_vs_wildtype_light_16h_unfiltered.csv") %>%
as_tibble() %>%
arrange(desc(log2FoldChange))
gene_list.TLRdarkWT16h <- res.TLRdarkWT16h$log2FoldChange
names(gene_list.TLRdarkWT16h) <- res.TLRdarkWT16h$ENSG
set.seed(4735)
gseGO.TLRdarkWT16h <- gseGO(gene_list.TLRdarkWT16h,
OrgDb = org.Hs.eg.db,
keyType = "ENSEMBL",
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = TRUE)## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...gseGO.TLRdarkWT16h <- setReadable(gseGO.TLRdarkWT16h, "org.Hs.eg.db", "ENSEMBL")
gseGO.TLRdarkWT16h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism Homo sapiens
## #...@setType BP
## #...@keytype ENSEMBL
## #...@geneList Named num [1:32287] 7.73 3.18 2.8 2.75 2.74 ...
## - attr(*, "names")= chr [1:32287] "ENSG00000174123" "ENSG00000162496" "ENSG00000139292" "ENSG00000225756" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...32 enriched terms found
## 'data.frame': 32 obs. of 11 variables:
## $ ID : chr "GO:0007059" "GO:0098813" "GO:0051983" "GO:0000819" ...
## $ Description : chr "chromosome segregation" "nuclear chromosome segregation" "regulation of chromosome segregation" "sister chromatid segregation" ...
## $ setSize : int 404 297 125 220 53 72 61 58 47 47 ...
## $ enrichmentScore: num 0.629 0.656 0.757 0.679 0.837 ...
## $ NES : num 1.77 1.81 1.89 1.8 1.9 ...
## $ pvalue : num 4.00e-06 8.32e-06 1.29e-05 1.86e-05 4.33e-05 ...
## $ p.adjust : num 0.0193 0.02 0.0208 0.0224 0.0232 ...
## $ qvalue : num 0.0184 0.0191 0.0198 0.0214 0.0221 ...
## $ rank : num 3353 4105 3299 3995 3256 ...
## $ leading_edge : chr "tags=44%, list=10%, signal=40%" "tags=49%, list=13%, signal=43%" "tags=53%, list=10%, signal=48%" "tags=51%, list=12%, signal=45%" ...
## $ core_enrichment: chr "TTN/FBXO5/CENPK/DPF3/AURKB/NEK7/BRIP1/RMDN1/GOLGA8A/PRICKLE1/CENPI/DSN1/CHMP3/CENPQ/KNTC1/SKA2/RMI2/CENPU/KIF23"| __truncated__ "TTN/FBXO5/CENPK/DPF3/AURKB/BRIP1/RMDN1/PRICKLE1/CENPI/DSN1/CHMP3/CENPQ/KNTC1/RMI2/KIF23/PSMC3IP/CDC6/PMF1/CDK1/"| __truncated__ "FBXO5/DPF3/AURKB/KNTC1/RMI2/CDC6/CDK1/RAD18/HASPIN/NCAPG/SPDL1/BUB1B/GEN1/SMC2/NCAPH/ZWINT/ANAPC5/TTK/NCAPG2/SP"| __truncated__ "TTN/FBXO5/CENPK/DPF3/AURKB/RMDN1/PRICKLE1/CENPI/DSN1/CHMP3/KNTC1/RMI2/KIF23/CDC6/CDK1/RACGAP1/HASPIN/NCAPG/SPDL"| __truncated__ ...
## #...Citation
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
## clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
## The Innovation. 2021, 2(3):100141write.csv(data.frame(gseGO.TLRdarkWT16h), file = paste0(out.dir, "gsea_go.csv"), row.names = F)
TLRdarkWT16h.go.plot <- nrow(gseGO.TLRdarkWT16h) > 0dotplot(gseGO.TLRdarkWT16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
cnetplot(gseGO.TLRdarkWT16h, showCategory = 5, color.params = list(foldChange = gene_list.TLRdarkWT16h),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
gseGO.TLRdarkWT16h.pwts <- pairwise_termsim(gseGO.TLRdarkWT16h)
emapplot(gseGO.TLRdarkWT16h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
KEGG
ids <- bitr(names(gene_list.TLRdarkWT16h), fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)## 'select()' returned 1:many mapping between keys and columnsids <- ids[!duplicated(ids$ENSEMBL),]
gene_list.TLRdarkWT16h.entrez <- gene_list.TLRdarkWT16h[names(gene_list.TLRdarkWT16h) %in% ids$ENSEMBL]
names(gene_list.TLRdarkWT16h.entrez) <- ids$ENTREZID
gene_list.TLRdarkWT16h.entrez <- gene_list.TLRdarkWT16h.entrez[!duplicated(names(gene_list.TLRdarkWT16h.entrez))]set.seed(5923)
gseKEGG.TLRdarkWT16h <- gseKEGG(gene_list.TLRdarkWT16h.entrez,
organism = "hsa",
keyType = "ncbi-geneid",
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = T)## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...gseKEGG.TLRdarkWT16h <- setReadable(gseKEGG.TLRdarkWT16h, "org.Hs.eg.db", "ENTREZID")
gseKEGG.TLRdarkWT16h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism hsa
## #...@setType KEGG
## #...@keytype ENTREZID
## #...@geneList Named num [1:22414] 7.73 3.18 2.8 2.75 2.74 ...
## - attr(*, "names")= chr [1:22414] "81793" "9249" "8549" "138948" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...2 enriched terms found
## 'data.frame': 2 obs. of 11 variables:
## $ ID : chr "hsa05204" "hsa00190"
## $ Description : chr "Chemical carcinogenesis - DNA adducts" "Oxidative phosphorylation"
## $ setSize : int 38 131
## $ enrichmentScore: num -0.789 0.641
## $ NES : num -1.89 1.74
## $ pvalue : num 0.000103 0.000285
## $ p.adjust : num 0.0346 0.0476
## $ qvalue : num 0.0316 0.0435
## $ rank : num 1887 3982
## $ leading_edge : chr "tags=32%, list=8%, signal=29%" "tags=56%, list=18%, signal=46%"
## $ core_enrichment: chr "GSTM4/CBR1/CYP1A1/GSTA4/GSTM3/GSTM2/MGST1/AKR1C2/EPHX1/CYP1B1/HSD11B1/PTGS2" "ATP6V0E2/ATP5F1D/COX3/COX2/ATP5F1C/NDUFB8/ATP5PO/NDUFA10/ATP5MG/NDUFB1/ATP5F1B/NDUFA11/ATP5F1A/COX6C/COX7B/NDUF"| __truncated__
## #...Citation
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
## clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
## The Innovation. 2021, 2(3):100141write.csv(data.frame(gseKEGG.TLRdarkWT16h), file = paste0(out.dir, "gsea_kegg.csv"), row.names = F)
TLRdarkWT16h.kegg.plot <- nrow(gseKEGG.TLRdarkWT16h) > 0dotplot(gseKEGG.TLRdarkWT16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
cnetplot(gseKEGG.TLRdarkWT16h, showCategory = 5, color.params = list(foldChange = gene_list.TLRdarkWT16h.entrez),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
gseKEGG.TLRdarkWT16h.pwts <- pairwise_termsim(gseKEGG.TLRdarkWT16h)
emapplot(gseKEGG.TLRdarkWT16h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
pathview(gene.data = gene_list.TLRdarkWT16h.entrez, pathway.id = gseKEGG.TLRdarkWT16h@result$ID[1], species = "hsa")## 'select()' returned 1:1 mapping between keys and columns## Info: Working in directory C:/Users/Daniel Zimmermann/Desktop/AnS_DESeq## Info: Writing image file hsa05204.pathview.pngfilename <- paste0(gseKEGG.TLRdarkWT16h@result$ID[1], ".pathview.png")
file.rename(filename, paste0(out.dir, filename))## [1] TRUEknitr::include_graphics(paste0(out.dir, filename))
Disease ontology
set.seed(6246)
gseDO.TLRdarkWT16h <- gseDO(gene_list.TLRdarkWT16h.entrez,
eps = 0,
maxGSSize = 500,
minGSSize = 15,
seed = T)## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...gseDO.TLRdarkWT16h <- setReadable(gseDO.TLRdarkWT16h, "org.Hs.eg.db", "ENTREZID")
gseDO.TLRdarkWT16h## #
## # Gene Set Enrichment Analysis
## #
## #...@organism Homo sapiens
## #...@setType DO
## #...@keytype ENTREZID
## #...@geneList Named num [1:22414] 7.73 3.18 2.8 2.75 2.74 ...
## - attr(*, "names")= chr [1:22414] "81793" "9249" "8549" "138948" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...3 enriched terms found
## 'data.frame': 3 obs. of 11 variables:
## $ ID : chr "DOID:9974" "DOID:2559" "DOID:974"
## $ Description : chr "drug dependence" "opiate dependence" "upper respiratory tract disease"
## $ setSize : int 18 17 154
## $ enrichmentScore: num 0.941 0.941 0.651
## $ NES : num 1.89 1.86 1.8
## $ pvalue : num 5.02e-06 1.57e-05 1.46e-04
## $ p.adjust : num 0.00482 0.00754 0.04667
## $ qvalue : num 0.00445 0.00697 0.04313
## $ rank : int 21 21 663
## $ leading_edge : chr "tags=11%, list=0%, signal=11%" "tags=12%, list=0%, signal=12%" "tags=6%, list=3%, signal=6%"
## $ core_enrichment: chr "KCNJ6/ABAT" "KCNJ6/ABAT" "TLR10/KCNJ6/BDNF/EDN1/CD40/MCAM/NLRP3/TLR4/TGFBR2/ADM"
## #...Citation
## Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an
## R/Bioconductor package for Disease Ontology Semantic and Enrichment
## analysis. Bioinformatics 2015, 31(4):608-609write.csv(data.frame(gseDO.TLRdarkWT16h), file = paste0(out.dir, "gsea_do.csv"), row.names = F)
TLRdarkWT16h.DO.plot <- nrow(gseDO.TLRdarkWT16h) > 0dotplot(gseDO.TLRdarkWT16h, showCategory = 15, size = NULL, split = ".sign", font.size = 9, label_format = 40) +
facet_grid(~.sign) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
cnetplot(gseDO.TLRdarkWT16h, showCategory = 5, color.params = list(foldChange = gene_list.TLRdarkWT16h.entrez),
cex.params = list(category_label = 0.6, gene_label = 0.4)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", limits = c(-2, 2), name = "Fold Change", oob = scales::squish)## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
gseDO.TLRdarkWT16h.pwts <- pairwise_termsim(gseDO.TLRdarkWT16h)
emapplot(gseDO.TLRdarkWT16h.pwts, force = 1.5,
cex.params = list(category_label = 0.8)) +
scale_fill_gradient(low = "red", high = "white", limits = c(0, 0.05), name = "adj. p-value")## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.